In [4]:
import cmocean as cmo
import numpy as np
import scipy.constants as spc
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation
from IPython.display import display, Math, HTML
from matplotlib.animation import FuncAnimation
plt.rcParams.update(
    {
        'mathtext.fontset':'cm', 
        'font.family':'serif', 
        'font.sans-serif':'Times New Roman', 
        'figure.figsize':[9,4], 
        'figure.titlesize':22,
        "figure.dpi":90, 
        'savefig.dpi':90,
        'axes.titlesize':20, 
        'axes.labelsize':18, 
        'axes.titley': 1.0, 
        'axes.titlepad': 5.0,
        'axes.edgecolor':'black', 
        'axes.grid': False, 
        'grid.alpha': .5,
        'xtick.labelsize':14,'ytick.labelsize':14,
        'xtick.major.size':6,'ytick.major.size':6,
        'xtick.major.width':1.25, 'ytick.major.width':1.25, 
        'xtick.direction':'inout','ytick.direction':'inout',
        'xtick.top':False, 'ytick.right':False,
        'legend.title_fontsize':14, 'legend.fontsize':14,
        'legend.borderaxespad': 1, 'legend.borderpad': 0.5,
        'legend.framealpha': 1,
        'legend.handleheight': 0.5, 'legend.handlelength': 2.0, 'legend.handletextpad':0.5,
        'legend.labelspacing': 0.25,
        'legend.fancybox':False,
        'legend.edgecolor': '0',
        'legend.frameon': True,
        'legend.markerscale': 1.25,
        'animation.embed_limit':2**128,
        'animation.html': 'jshtml'
    }
)
pi = spc.pi
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
default_plots()

def periodic(f):
    f[0,:]  = f[-2,:]
    f[-1,:] = f[2,:]
    f[:,0]  = f[:,-2]
    f[:,-1] = f[:,2]
    return

def fixed(f):
    f[0,:]  = 0
    f[-1,:] = 0
    f[:,0]  = 0
    f[:,-1] = 0
    return

c = 1
x = np.linspace(0,1,128)
y = np.linspace(0,1,128)
xx, yy = np.meshgrid(x,y)

u0 = np.zeros(y.shape+x.shape)

steps = 500

dx = x[1]-x[0]
dy = y[1]-y[0]
dt = (dx/c)/2

t = np.linspace(0,steps*dt,steps+1)
u = np.zeros(t.shape + u0.shape)
u_t = np.zeros(t.shape + u0.shape)

u[0] = u0/10
u_t[0] = np.zeros_like(u0)

rx = (c/dx)**2
ry = (c/dy)**2

# def source(f, n_t):
#     return np.cos(f*pi*n_t/steps)**2


xL = .5-dx
xR = .5+dx
yB = .5-3*dy
yT = .5+3*dy
mask = ~(((xx < xR)*(xx > xL)*(yy>yT)) + ((xx < xR)*(xx > xL)*(yy<yB)))
source = (xx-.4)**2 + (yy-.5)**2 < .002

boundary = 'absorbing'
src_int = 40
stp_int = 2

if boundary == 'periodic':
    for n in range(steps):
        if (n%src_int==0)*(n < steps//stp_int):
            u[n] += source
            
        u = u*mask
        
        periodic(u[n])
        
        u_t[n+1,1:-1,1:-1] = u_t[n,1:-1,1:-1] \
                    + dt*rx * (u[n,2:,1:-1] + u[n,:-2,1:-1] - 2*u[n,1:-1,1:-1]) \
                    + dt*ry * (u[n,1:-1,2:] + u[n,1:-1,:-2] - 2*u[n,1:-1,1:-1])
        u[n+1,1:-1,1:-1] = u[n,1:-1,1:-1] + dt*u_t[n+1,1:-1,1:-1]
        
if boundary == 'fixed':
    for n in range(steps):
        if (n%src_int==0)*(n < steps//stp_int):
            u[n] += source
            
        u = u*mask
        
        fixed(u[n])
        
        u_t[n+1,1:-1,1:-1] = u_t[n,1:-1,1:-1] \
                    + dt*rx * (u[n,2:,1:-1] + u[n,:-2,1:-1] - 2*u[n,1:-1,1:-1]) \
                    + dt*ry * (u[n,1:-1,2:] + u[n,1:-1,:-2] - 2*u[n,1:-1,1:-1])
        u[n+1,1:-1,1:-1] = u[n,1:-1,1:-1] + dt*u_t[n+1,1:-1,1:-1]
        
if boundary == 'absorbing':
    for n in range(steps):
        if (n%src_int==0)*(n < steps//stp_int):
            u[n] += source
        u = u*mask
        
        u[n+1,1:-1,0] = (1-c*dt/dx)*u[n,1:-1,0] + (c*dt/dx)*u[n,1:-1,1]
        u[n+1,1:-1,-1] = (1-c*dt/dx)*u[n,1:-1,-1] + (c*dt/dx)*u[n,1:-1,-2]
        u[n+1,0,1:-1] = (1-c*dt/dy)*u[n,0,1:-1] + (c*dt/dy)*u[n,1,1:-1]
        u[n+1,-1,1:-1] = (1-c*dt/dy)*u[n,-1,1:-1] + (c*dt/dy)*u[n,-2,1:-1]
        
        u_t[n+1,1:-1,1:-1] = u_t[n,1:-1,1:-1] \
                    + dt*rx * (u[n,2:,1:-1] + u[n,:-2,1:-1] - 2*u[n,1:-1,1:-1]) \
                    + dt*ry * (u[n,1:-1,2:] + u[n,1:-1,:-2] - 2*u[n,1:-1,1:-1])
        u[n+1,1:-1,1:-1] = u[n,1:-1,1:-1] + dt*u_t[n+1,1:-1,1:-1]
        
def animate2D_mesh(u, t, fps = 24, frames = 96, cmap = cmo.cm.deep):
    t_steps = int(len(t)-1)
    i_frames = np.arange(0, t_steps, int(t_steps/frames))
    fig, ax = plt.subplots(figsize = (7,7), constrained_layout = True)
    image = ax.imshow(u[0], cmap = cmap, vmin = np.min(u[0])/2, vmax = np.max(u[0])/2)
    ax.axis(False)
    def animate(i):
        image.set_array(u[i])
    plt.close()
    anim = FuncAnimation(fig, animate, interval = int(1e3/fps), frames = i_frames)
    return anim
In [5]:
fps = 24
frames = 144
anim = animate2D_mesh(u, t, fps = fps, frames = frames)
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
display(anim)